home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_10_05 / 1005036a < prev    next >
Encoding:
Text File  |  1992-03-21  |  3.1 KB  |  78 lines
  1. #include <math.h>
  2.  
  3. /* Adaptive quadrature routine.
  4.    Called with:
  5.         lft:            left endpoint
  6.         rt:             right endpoint
  7.         f:              pointer to function to integrate
  8.         err:            pointer to allowable error
  9.         p1:             function value at left endpoint
  10.         p2:             function value at center of interval
  11.         p3:             function value at right endpoint
  12.  
  13.    Returns
  14.         return value    evaluated integral
  15.         err             pointer to upper limit of actual error
  16. */
  17.  
  18. static double adaptive2 (double lft, double rt, double (*f)(double x), double *err,
  19.                   double p1, double p3, double p5)
  20. {
  21.    double p2, p4;                      /* Function values at pts 2 & 4. */
  22.    double h;                           /* Width of this interval. */
  23.    double error;                       /* Actual error on this interval. */
  24.    double err1, err2;                  /* Errors of the sub-intervals. */
  25.    double int3,                        /* Value of the integral with Simpson */
  26.           int5;                        /* 3- and 5-point approximations. */
  27.  
  28.    h = rt - lft;                       /* Get the width of the interval. */
  29.                               /* Get the 2 approximations to the integral. */
  30.    p2 = (*f)(lft + h/4.0);
  31.    p4 = (*f)(rt - h/4.0);
  32.  
  33.    int3 = h * (p1 + 4.0 * p3 + p5) / 6.0;
  34.    int5 = h * (p1 + 4.0 * p2 + 2.0 * p3 + 4.0 * p4 + p5) / 12.0;
  35.  
  36.    error = fabs (int3 - int5) / 15.0;     /* Find the actual error. */
  37.    if (error < *err)                      /* If within tolerance, */
  38.    {
  39.       *err = error;                       /* copy the actual error to caller, */
  40.       return (int5);                      /* return the more accurate approx. */
  41.    }
  42.    else                                   /* Otherwise, */
  43.    {
  44.       err1 = err2 = *err / 2.0;           /* integrate the two half-regions, */
  45.       int5 = adaptive2 (lft, lft + h/2.0, f, &err1, p1, p2, p3);
  46.       int5 += adaptive2 (lft + h/2.0, rt, f, &err2, p3, p4, p5);
  47.       *err = err1 + err2;                 /* sum the real errors, */
  48.       return (int5);                      /* and return the integral. */
  49.    }
  50. }
  51.  
  52. /* Adaptive quadrature interface.
  53.    Called with:
  54.         lft:            left endpoint
  55.         rt:             right endpoint
  56.         f:              pointer to function to integrate
  57.         err:            pointer to allowable error
  58.  
  59.    Returns
  60.         return value    evaluated integral
  61.         err             pointer to upper limit of actual error
  62. */
  63. double adaptive (double lft, double rt, double (*f)(double x), double *err)
  64. {
  65.    double p1, p2, p3;                  /* Initial function evaluations. */
  66.    double h;                           /* Width of this interval. */
  67.  
  68.    h = rt - lft;                       /* Get the width of the interval. */
  69.  
  70.    p1 = (*f)(lft);                     /* Pre-compute the initial function */
  71.    p2 = (*f)(lft + h/2.0);             /* values. */
  72.    p3 = (*f)(rt);
  73.  
  74.                                      /* Call adaptive2 to do the real work. */
  75.    return (adaptive2 (lft, rt, f, err, p1, p2, p3));
  76. }
  77.  
  78.